Observed log-chlorophyll at representative station for the St. Lucie Estuary
library(tidyverse)
library(lubridate)
library(mgcv)
library(plotly)
library(WRTDStidal)
library(gridExtra)
source('R/funcs.R')
# format the data to model
data(sl_dat)
sl_mod <- sl_dat %>%
rename(date = Date) %>%
mutate(
doy = yday(date),
dec_time = decimal_date(date),
yr = year(date),
mo = month(date, label = T)
) %>%
mutate(
flo = sal,
lnchl = log(1 + chl)
) %>%
select(-sal)
# plot, all
p <- ggplot(sl_mod, aes(x = date, y = lnchl)) +
geom_line() +
theme_bw()
ggplotly(p)
# boxplot, by years
p <- ggplot(sl_mod, aes(x = yr, y = lnchl)) +
geom_boxplot() +
theme_bw()
ggplotly(p)
# boxplot, by month
p <- ggplot(sl_mod, aes(x = mo, y = lnchl)) +
geom_boxplot() +
theme_bw()
ggplotly(p)
Some simple GAMs to explore annual, seasonal trends.
# smooths to evaluate
smths <- c(
"s(dec_time, bs = 'tp')",
"s(doy, bs = 'cc')",
"te(dec_time, doy, bs = c('tp', 'cc'))"
)
# get all combinations of smoothers to model, one to many
frms <- list()
for(i in seq_along(smths)){
frm <- combn(smths, i) %>%
apply(2, function(x){
paste(x, collapse = ' + ') %>%
paste('lnchl ~ ', .) %>%
formula
})
frms <- c(frms, frm)
}
# create models from smooth formula combinations
mods <- map(frms, function(frm){
gam(frm,
knots = list(doy = c(1, 366)),
data = sl_mod,
na.action = na.exclude
)
})
names(mods) <- paste0('mod', seq_along(mods))
Summary stats of annual, seasonal models:
# smoother stats of GAMs
map(mods, ~ summary(.x)$s.table %>% data.frame %>% rownames_to_column('smoother')) %>%
enframe %>%
unnest %>%
kable
| name | smoother | edf | Ref.df | F | p.value |
|---|---|---|---|---|---|
| mod1 | s(dec_time) | 1.6461621 | 2.048370 | 9.7038674 | 0.0000678 |
| mod2 | s(doy) | 4.0986733 | 8.000000 | 6.8154674 | 0.0000000 |
| mod3 | te(dec_time,doy) | 9.9947633 | 12.933682 | 5.5626179 | 0.0000000 |
| mod4 | s(dec_time) | 1.9430775 | 2.429282 | 8.9065355 | 0.0000567 |
| mod4 | s(doy) | 4.1145399 | 8.000000 | 7.0520936 | 0.0000000 |
| mod5 | s(dec_time) | 5.9855784 | 7.226382 | 2.1270530 | 0.0369088 |
| mod5 | te(dec_time,doy) | 12.3369743 | 15.000000 | 4.0303044 | 0.0000000 |
| mod6 | s(doy) | 4.0464881 | 8.000000 | 6.3773815 | 0.0000000 |
| mod6 | te(dec_time,doy) | 2.4100580 | 3.304212 | 6.7570411 | 0.0001249 |
| mod7 | s(dec_time) | 1.9175106 | 2.397382 | 9.0127722 | 0.0000550 |
| mod7 | s(doy) | 4.1236170 | 8.000000 | 7.0552277 | 0.0000000 |
| mod7 | te(dec_time,doy) | 0.0000065 | 15.000000 | 0.0000004 | 0.3651312 |
# summary stats of GAMs
map(mods, ~ data.frame(
AIC = AIC(.x),
R2 = summary(.x)$r.sq)) %>%
enframe %>%
unnest %>%
kable
| name | AIC | R2 |
|---|---|---|
| mod1 | 683.2360 | 0.0500261 |
| mod2 | 653.5304 | 0.1302894 |
| mod3 | 643.6646 | 0.1668405 |
| mod4 | 634.5680 | 0.1787925 |
| mod5 | 644.2521 | 0.1835579 |
| mod6 | 634.4836 | 0.1798639 |
| mod7 | 634.5693 | 0.1787532 |
pred_dat <- sl_mod
# predictions
sl_res <- map(mods, function(x){
pred_dat %>%
mutate(
pred = predict(x, newdata = pred_dat)
)
}) %>%
enframe('mods') %>%
unnest
# plot
p <- ggplot(sl_res, aes(x = date)) +
geom_point(data = sl_mod, aes(y = lnchl), size = 0.5) +
geom_line(aes(y = pred, colour = mods)) +
theme_bw() +
theme(
legend.position = 'top',
legend.title = element_blank()
)
ggplotly(p)
# plot
p <- ggplot(sl_res, aes(x = doy, group = factor(yr), colour = yr)) +
geom_line(aes(y = pred)) +
theme_bw() +
theme(
legend.position = 'top',
legend.title = element_blank()
) +
facet_wrap(~ mods, ncol = 2)
ggplotly(p)
Adding flow data to the model:
# smooths to evaluate
smths <- c(
"s(dec_time, bs = 'tp')",
"s(doy, bs = 'cc')",
"s(flo, bs = 'tp')",
"te(flo, doy, bs = c('tp', 'cc'))",
"te(flo, dec_time, bs = c('tp', 'tp'))",
"te(dec_time, doy, bs = c('tp', 'cc'))",
"te(dec_time, doy, flo, bs = c('tp', 'cc', 'tp'))"
)
# get all combinations of smoothers to model, one to many
frms <- list()
for(i in seq_along(smths)){
frm <- combn(smths, i) %>%
apply(2, function(x){
paste(x, collapse = ' + ') %>%
paste('lnchl ~ ', .) %>%
formula
})
frms <- c(frms, frm)
}
# create models from smooth formula combinations
mods2 <- map(frms, function(frm){
gam(frm,
knots = list(doy = c(1, 366)),
data = sl_mod,
na.action = na.exclude
)
})
names(mods2) <- paste0('mod', seq_along(mods2))
Summary stats of best year/season model, year/season/flow model
# best model with only season, year
best1 <- map(mods, AIC) %>%
unlist %>%
which.min %>%
mods[[.]]
# best model with season, year, flow
best2 <- map(mods2, AIC) %>%
unlist %>%
which.min %>%
mods2[[.]]
best <- list(best1 = best1, best2 = best2)
# smoother stats of GAMs
map(best, ~ summary(.x)$s.table %>% data.frame %>% rownames_to_column('smoother')) %>%
enframe %>%
unnest %>%
kable
| name | smoother | edf | Ref.df | F | p.value |
|---|---|---|---|---|---|
| best1 | s(doy) | 4.0464881 | 8.000000 | 6.3773815 | 0.0000000 |
| best1 | te(dec_time,doy) | 2.4100580 | 3.304212 | 6.7570411 | 0.0001249 |
| best2 | s(doy) | 3.2726969 | 8.000000 | 5.4646634 | 0.0000000 |
| best2 | s(flo) | 2.1261492 | 2.532475 | 1.1775041 | 0.2361336 |
| best2 | te(flo,doy) | 0.9077802 | 13.000000 | 0.1933511 | 0.0154606 |
| best2 | te(flo,dec_time) | 7.4271874 | 20.000000 | 1.0880301 | 0.0000178 |
| best2 | te(dec_time,doy) | 0.5274598 | 8.000000 | 0.0673205 | 0.2037872 |
| best2 | te(dec_time,doy,flo) | 12.3847743 | 48.000000 | 0.6027438 | 0.0000674 |
# summary stats of GAMs
map(best, ~ data.frame(
AIC = AIC(.x),
R2 = summary(.x)$r.sq)) %>%
enframe %>%
unnest %>%
kable
| name | AIC | R2 |
|---|---|---|
| best1 | 634.4836 | 0.1798639 |
| best2 | 502.3692 | 0.3272847 |
# predictions
sl_res2 <- map(best, function(x){
pred_dat %>%
mutate(
pred = predict(x, newdata = pred_dat)
)
}) %>%
enframe('mods') %>%
unnest
# plot
p <- ggplot(sl_res2, aes(x = date)) +
geom_point(data = sl_mod, aes(y = lnchl), size = 0.5) +
geom_line(aes(y = pred, colour = mods)) +
theme_bw() +
theme(
legend.position = 'top',
legend.title = element_blank()
)
ggplotly(p)
ptheme <- theme(
axis.title.x = element_blank(),
axis.title.y = element_blank()
)
cols <- 'Spectral'
pb1 <- dynagam(best1, pred_dat, ncol = 1, col_vec = cols) +
ptheme +
theme(legend.position = 'none') +
ggtitle('Best 1')
pb2 <- dynagam(best2, pred_dat, ncol = 1, col_vec = cols) +
ptheme +
ggtitle('Best2')
pleg <- g_legend(pb2)
pb2 <- pb2 +
theme(legend.position = 'none')
grid.arrange(
pleg,
arrangeGrob(pb1, pb2, ncol = 2, bottom = 'lnQ', left = 'lnchl'),
ncol = 1,
heights = c(0.1, 1)
)